﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Hydra
{
  class Hydra
  {
    public int n, m;
    public double[] u1, u2;
    public double lambda1, lambda2, ro, mu, nu, k, c;
    public double maxT, maxX;
    public double h, dt;
    public double u1w1, u2w2, u1w2, u2w1;
    public double p;

    public void hydraConstr(int N, int M,double MaxT,double MaxX, double Lambda1, double Lambda2, double Ro, double Mu, double Nu, double K, double C, double P)
    {
      maxT = MaxT;
      maxX = MaxX;
      h = maxX / N;
      dt = maxT / M;
      p = P;
      lambda1 = Lambda1;
      lambda2 = Lambda2;
      ro = Ro;
      mu = Mu;
      nu = Nu;
      k = K;
      c = C;
      n = N;
      m = M;
      u1 = new double[(n + 1) * (m + 1)];
      u2 = new double[(n + 1) * (m + 1)];
    }
    // 0 <= i <= n
    // 0 <= j <= m
    public int getIdx(int i, int j, int n)
    {
      return (j * (n + 1) + i);
    }
    public void hydraSetFreq(double U1w1, double U1w2, double U2w1, double U2w2)
    {
      u1w1 = U1w1;
      u1w2 = U1w2;
      u2w1 = U2w1;
      u2w2 = U2w2;
    }
    public void useBoundarys()
    {
      // Нижняя граница
      for (int i = 0; i <= n; i++)
      {
        int j = 0;
        int idx = getIdx(i, j, n);
        u1[idx] = Math.Cos(u1w1 * i * h - u1w2 * j * dt);
        u2[idx] = Math.Cos(u2w1 * i * h - u2w2 * j * dt);
      }

      // Левая граница
      for (int j = 0; j <= m; j++)
      {
        int i = 0;
        int idx = getIdx(i, j, n);
        u1[idx] = Math.Cos(u1w1 * i * h - u1w2 * j * dt);
        u2[idx] = Math.Cos(u1w1 * i * h - u1w2 * j * dt);
      }

      // Правая граниа
      for (int j = 0; j <= m; j++)
      {
        int i = n;
        int idx = getIdx(i, j, n);
        u1[idx] = Math.Cos(u1w1 * i * h - u1w2 * j * dt);
        u2[idx] = Math.Cos(u1w1 * i * h - u1w2 * j * dt);
      }
    }    
    // Методы решения.
    // явный
    public void solveExplicit2gu()
    {
      u1[getIdx(1, 0, n)] = u1[getIdx(0, 0, n)];
      u1[getIdx(n - 1, 0, n)] = u1[getIdx(n, 0, n)];
      u2[getIdx(1, 0, n)] = u2[getIdx(0, 0, n)];
      u2[getIdx(n - 1, 0, n)] = u2[getIdx(n, 0, n)];
      for (int j = 0; j < m; j++)
      {
        for (int i = 1; i < n; i++)
        {

          u1[getIdx(i, j + 1, n)] = ((u1[getIdx(i - 1, j, n)] - 2 * u1[getIdx(i, j, n)] + u1[getIdx(i + 1, j, n)]) * lambda1 / (h * h)
            + ro - mu * u1[getIdx(i, j, n)] + k * u1[getIdx(i, j, n)] * u1[getIdx(i, j, n)] / u2[getIdx(i, j, n)]) * dt + u1[getIdx(i, j, n)];
          u2[getIdx(i, j + 1, n)] = ((u2[getIdx(i - 1, j, n)] - 2 * u2[getIdx(i, j, n)] + u2[getIdx(i + 1, j, n)]) * lambda2 / (h * h)
            - nu * u2[getIdx(i, j, n)] + c * u1[getIdx(i, j, n)] * u1[getIdx(i, j, n)]) * dt + u2[getIdx(i, j, n)];
        }
        u1[getIdx(0, j + 1, n)] = u1[getIdx(1, j + 1, n)];
        u1[getIdx(n, j + 1, n)] = u1[getIdx(n - 1, j + 1, n)];
        u2[getIdx(0, j + 1, n)] = u2[getIdx(1, j + 1, n)];
        u2[getIdx(n, j + 1, n)] = u2[getIdx(n - 1, j + 1, n)];
      }
      //u1[getIdx(1, m, n)] = u1[getIdx(0, m, n)];
      //u1[getIdx(n - 1, m, n)] = u1[getIdx(n, m, n)];
      //u2[getIdx(1, m, n)] = u2[getIdx(0, m, n)];
      //u2[getIdx(n - 1, m, n)] = u2[getIdx(n, m, n)];
    }
    // неявный
    public void solveImplicit2gu()
    {
      // Учитываем 2 тип ГУ
      for (int j = 0; j <= m; j++)
      {
        u1[getIdx(1, j, n)] = u1[getIdx(0, j, n)];
        u1[getIdx(n - 1, j, n)] = u1[getIdx(n, j, n)];

        u2[getIdx(1, j, n)] = u2[getIdx(0, j, n)];
        u2[getIdx(n - 1, j, n)] = u2[getIdx(n, j, n)];
      }


      // Решаем СЛАУ, с 3-х диагональной матрицей методом прогонки.
      double[] Answer = new double[n + 1];
      double[] alpha = new double[n + 1];
      double[] beta = new double[n + 1];
      double[] A = new double[n];
      double[] B = new double[n];
      double[] C = new double[n];
      double[] Phi = new double[n];
      // цикл по слоям
      for (int j = 1; j <= m; j++)
      {
        #region обновление u1
        // Заполняем коэффициенты для прогонки.
        for (int i = 1; i < n; i++)
        {
          C[i] = -(1.0 / dt + 2 * lambda1 / (h * h) + mu);
          Phi[i] = -(u1[getIdx(i, j - 1, n)] / dt + ro + k 
            * u1[getIdx(i, j - 1, n)] * u1[getIdx(i, j - 1, n)] 
            / u2[getIdx(i, j - 1, n)]);
          A[i] = -lambda1 / (h * h);
          B[i] = -lambda1 / (h * h);
        }

        Progonka(1, 1, 0, 0, alpha, beta, A, B, C, Phi, Answer);

        for (int i = 0; i <= n; i++)
        {
          u1[getIdx(i, j, n)] = Answer[i];
        }
        #endregion 

        #region обновление u2
        // Заполняем коэффициенты для прогонки.
        for (int i = 1; i < n; i++)
        {
          C[i] = -(1.0 / dt + 2 * lambda2 / (h * h) + nu);
          Phi[i] = -(u2[getIdx(i, j - 1, n)] / dt + c 
            * u1[getIdx(i, j - 1, n)] * u1[getIdx(i, j - 1, n)]);
          A[i] = -lambda2 / (h * h);
          B[i] = -lambda2 / (h * h);
        }

        Progonka(1, 1, 0, 0, alpha, beta, A, B, C, Phi, Answer);

        for (int i = 0; i <= n; i++)
        {
          u2[getIdx(i, j, n)] = Answer[i];
        }
        #endregion 
      }
    }
    // схема с весом.
    public void solveWeight()
    {
      // Учитываем 2 тип ГУ
      u1[getIdx(1, 0, n)] = u1[getIdx(0, 0, n)];
      u1[getIdx(n - 1, 0, n)] = u1[getIdx(n, 0, n)];
      u2[getIdx(1, 0, n)] = u2[getIdx(0, 0, n)];
      u2[getIdx(n - 1, 0, n)] = u2[getIdx(n, 0, n)];

      // Решаем СЛАУ, с 3-х диагональной матрицей методом прогонки.
      double[] Answer = new double[n + 1];
      double[] alpha = new double[n + 1];
      double[] beta = new double[n + 1];
      double[] A = new double[n];
      double[] B = new double[n];
      double[] C = new double[n];
      double[] Phi = new double[n];

      // цикл по слоям
      for (int j = 1; j <= m; j++)
      {
        #region обновление u1
        // Заполняем коэффициенты для прогонки.
        for (int i = 1; i < n; i++)
        {
          C[i] = -(1.0 / dt + 2 * lambda1 * (1 - p) / (h * h) + mu);
          Phi[i] = -(u1[getIdx(i, j - 1, n)] / dt + ro + k
            * u1[getIdx(i, j - 1, n)] * u1[getIdx(i, j - 1, n)]
            / u2[getIdx(i, j - 1, n)] + lambda1 * p * (u1[getIdx(i - 1, j - 1, n)]
            - 2 * u1[getIdx(i, j - 1, n)] + u1[getIdx(i + 1, j - 1, n)]) 
            / (h * h));
          A[i] = -lambda1 * (1 - p) / (h * h);
          B[i] = -lambda1 * (1 - p) / (h * h);
        }

        Progonka(1, 1, 0, 0, alpha, beta, A, B, C, Phi, Answer);

        for (int i = 0; i <= n; i++)
        {
          u1[getIdx(i, j, n)] = Answer[i];
        }
        #endregion

        #region обновление u2
        // Заполняем коэффициенты для прогонки.
        for (int i = 1; i < n; i++)
        {
          C[i] = -(1.0 / dt + 2 * lambda2 * (1 - p) / (h * h) + nu);
          Phi[i] = -(u2[getIdx(i, j - 1, n)] / dt + c
            * u1[getIdx(i, j - 1, n)] * u1[getIdx(i, j - 1, n)])
            + lambda2 * p * (u2[getIdx(i - 1, j - 1, n)]
            - 2 * u2[getIdx(i, j - 1, n)] + u2[getIdx(i + 1, j - 1, n)]) 
            / (h * h);
          A[i] = -lambda2 * (1 - p) / (h * h);
          B[i] = -lambda2 * (1 - p) / (h * h);
        }

        Progonka(1, 1, 0, 0, alpha, beta, A, B, C, Phi, Answer);

        for (int i = 0; i <= n; i++)
        {
          u2[getIdx(i, j, n)] = Answer[i];
        }
        #endregion
      }
    }

    // Метод прогонки, для 3-х диагональной матрицы, обозначения из лекций
    // Стронгиной Н. Р.
    void Progonka(double kappa1, double kappa2, double mu1, double mu2,
      double[] alpha, double[] beta, double[] A, double[] B, double[] C,
      double[] Phi, double[] Answer)
    {
      alpha[1] = kappa1;
      beta[1] = mu1;
      for (int i = 1; i < n; i++)
      {
        alpha[i + 1] = B[i] / (C[i] - alpha[i] * A[i]);
        beta[i + 1] = (Phi[i] + A[i] * beta[i]) / (C[i] - alpha[i] * A[i]);
      }
      Answer[n] = (mu2 + kappa2 * beta[n]) / (1 - kappa2 * alpha[n]);
      for (int i = n - 1; i >= 0; i--)
      {
        Answer[i] = alpha[i + 1] * Answer[i + 1] + beta[i + 1];
      }
    }
  }
}
